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We present a novel approach to determine the onset of contact between a tip and a surface. The 
real contact area depending on the distance is calculated using Bader’s quantum theory of atoms in 
molecules. The jump to contact, which is often observed in atomic force microscopy experiments, is 
used as an indicator for the initial point of contact, which in turn is defined by atomic relaxations and 
thus without the need of external parameters. Within our approach the contact area is estimated by 
evaluating the zero flux surfaces between the touching Bader-atoms, where the necessary electronic 
density cutoff for the Bader-partitioning is calculated to depend on the initial point of contact. Our 
proposed approach is therefore completely ah-initio and we are able to define and calculate the real 
area of contact without imposing restrictions or free parameters. As a prototype system we choose 
a tip made of a ten atom tungsten pyramid above a moire layer of graphene on an fee iridium (111) 
substrate. We find that the contact area depends exponentially on the effective distance between 
the tip apex and the surface atom directly below within the atomically relaxed nanosystem. 

PACS numbers: 71.15.Mb, 68.37.Ef,68.37.Ps,73.63.Rt,68.55.ap 


INTRODUCTION 

Historically the introduction of the concept of a real 
contact area in 1939 by Bowden and Tabor, which is sub¬ 
stantially smaller than the nominal one, was a huge step 
forward in understanding the laws of friction [T]. Since 
then a multitude of methods have been used to predict 
the pressure dependent real contact area for rough sur¬ 
faces. 

In 1881 Hertz laid the foundations of contact mechan¬ 
ics by describing the junctions between non-adhesive, ho¬ 
mogeneous elastic solids of simple shapes [2]. Nearly 
a hundred years later Johnson, Kendall, and Roberts 
(JKR) included short-range adhesion forces inside the 
contact which lead to larger contact areas compared to 
Hertz’s model [3]. In contrast, Derjaguin, Muller and 
Toporov (DMT) assumed that the Hertzian contact area 
remains undeformed while long-range adhesion forces are 
acting outside the contact zone [4]. Tabor was able to 
show that both the JKR and the DMT model can be 
viewed as limiting cases of a more general model, with 
JKR suitable for soft materials with large adhesion and 
DMT for hard materials with low adhesion [5]. Ad¬ 
ditional work on this unification was done by Maugis 
using Dugdale potentials [6]. The analytical solution 
by Maugis, however, produces rather cumbersome equa¬ 
tions, which can be approximated with high accuracy by 
the generalized transition equation derived by Carpick, 
Ogletree and Salmeron [7]. Their result also describes 
the transition between the DMT and the JKR models, 
but with simpler expressions which can be applied more 
straightforwardly to experimental data and only differs 
from the Maugis-Dugdale model within the unstable low 
load region. Generally these theories agree that the con¬ 
tact area A between a single sphere and a flat surface is 
a sublinear function of the load L, see e.g. the original 


Hertz prediction of A{L) ^ 1? I ^. These continuum me¬ 
chanics models are undoubtedly very successful and play 
an important role in both theoretical and experimental 
work in tribology. However, on the atomic scale, as tested 
for example in AFM experiments, the size of the the con¬ 
tact approaches the size of the involved atoms and thus 
models based on continuum mechanics are hardly appli¬ 
cable [8]. The apparent success of the Maugis-Dugdale 
model, which is still widely used to interpret atomic scale 
AFM experiments, can be attributed to its flexibility pro¬ 
vided by three fitting parameters [9]. Hence, atomistic 
methods are of great value in assessing the validity of the 
results and provide a much needed tool to increase the 
understanding of the real contact area. 

In recent studies and discussions a need for charac¬ 
terizing the contact area on an atom-by-atom basis was 
expressed and various strategies to estimate the number 
of atoms in contact were proposed However, it 

is not trivial to decide at which distance two atoms are in 
contact or how big the resulting contact area should be. 
One possibility is to define a certain inter-atomic distance 
ao below which contact should be established, however, 
this just shifts the problem to find the correct distance 
ao- A common method is to identify ao, and thus the 
onset of contact, as the beginning of repulsion between 
the two observed atoms mun]. To this end, classical 
molecular dynamics (MD) studies using Lennard-Jones 
potentials are often employed, sometimes without the at¬ 
tractive part if the surfaces of interest are non-adhesive. 
In our view, this method is not ideal for describing the 
onset of contact on the nanoscale for the following two 
reasons: (i) if only repulsive interaction is a sign of con¬ 
tact, the atoms in a solid or molecule at OK are not in 
contact with each other; (ii) in AFM experiments one 
often observes a “jump to contact”, where either the tip, 
or some part of the surface below or both jump towards 


2 


each other because of strong attractive interactions. If 
now the tip support is lowered further, the distance be¬ 
tween the surface and the tip apex might get even smaller 
as the chemical binding becomes stronger. It seems un¬ 
reasonable to argue that the tip is not in contact with 
the surface after the jump, just because the interaction 
is still attractive. 

In a study from 2009, Mo, Turner, and Szlufarska ex¬ 
amined the contact area between hydrogenated amor¬ 
phous carbon tips of up to 30 nm radius and a flat hy¬ 
drogenated diamond surface m- For this large scale fi¬ 
nite temperature (300 K) MD study they used a reactive 
empirical bond-order (REBO) potential [15] to model 
the chemical forces and added an analytical switching 
function to include van der Waals (vdW) like long-range 
forces. The multi-asperity picture of nanoscale contact 
presented in their publication relies on the assumption 
that contact is established between the atoms that are 
interacting chemically through the REBO potential while 
a much larger part of the tip is attracted to the surface 
via the vdW forces. These ideas are further discussed in 
a follow-up publication by Mo and Szlufarska m- Now 
an atomic contact area Aat is attributed to every chemi¬ 
cally interacting atom of the tip leading to the total real 
contact area ^reai = ^at^at? with A^at being the num¬ 
ber of involved atoms. Due to atomic scale roughness 
in the amorphous tip this real contact area may be sig¬ 
nificantly smaller than the expected contact area of a 
smooth asperity, Aasp, which is defined as the envelope 
over the contact points. Concluding, Mo, Turner and, 
Szlufarska pointed out how important atomic corruga¬ 
tions are for an accurate estimation of the real contact 
area, and they underlined the need for accurate computa¬ 
tional approaches. Nevertheless, a few points remain to 
be analyzed further. In a realistic, continuous potential it 
might be hard to distinguish between long-range disper¬ 
sion forces and chemical forces, thus making it difficult 
to define the number of atoms in contact. Eurt her more, 
the inherent assumption that all contributions from sin¬ 
gle atomic contacts Aat are of equal size might not hold 
in systems with varying local load. 

In the following, we will introduce an ah-initio ap¬ 
proach to estimate the distance-dependent real contact 
area between a tip and a surface. To prove the feasibility 
of our approach, we choose a realistic and thus rather 
complex system, which has previously been used to ex¬ 
plain contrast inversion between constant current scan¬ 
ning tunneling microscope and constant frequency AEM 
images of graphene moire on metals EHIH]. We employ 
van der Waals corrected density functional calculations 
and use Bader-partitioning of the electronic charge den¬ 
sity to identify accurate atomic volumes and surfaces in 
different chemical surroundings. 


COMPUTATIONAL METHODS 

The simulation cell consists of 4 layers of a 9 x 9 
iridium(lll) substrate covered by a 10 x 10 graphene 
layer forming a moire structure with an average sepa¬ 
ration of 3.42 A and a corrugation of 0.35 A. The lat¬ 
tice constants obtained with the optB86b-vdW func¬ 
tional m [20] . air = 2.735 A and acr = 2.465 A, are 
close to the experimental values of 2.71 A and 2.46 A, re¬ 
spectively. The mismatch of the structure is very small 
at lOacr — 9air = 0.015 A [18]. The tungsten AEM tip 
was modeled as a ten atom pyramid with one atom at the 
apex, four in the next layer and five at the top. The con¬ 
tact site studied was an on-top position in a top-hep re¬ 
gion of the moire structure. This means that the tip apex 
atom is positioned directly vertical over a carbon atom in 
a region where each carbon atom is either directly over 
an iridium atom or over an iridium hep position. The 
simulation cell, containing 534 atoms, is shown in fig¬ 
ure 12 Relaxations were allowed for the graphene layer 
and the bottom five atoms of the tip, keeping the iridium 
substrate and the top layer of the tip rigid at their initial 
relaxed positions. Relaxations of the iridium substrate 
during movement of the tip have been neglected due to 
the small binding energy (^80 meV per carbon atom) of 
the mainly physisorbed graphene layer, which makes any 
effect of the iridium substrate on the relaxations of the 
graphene unlikely. The topmost layer of the tungsten 
pyramid, on the other hand, needs to be held fixed to 
control the distance between the tip and the graphene 
sheet. 

All calculations were performed within density func¬ 
tional theory (DET) employing the Vienna Ab-Initio 
Simulation Package VASP [HHH] using the Projec¬ 
tor Augmented-Wave (PAW) method [25] |26]. To in¬ 
clude van der Waals (vdW) forces, which are relevant 
in this system, the optB86b-vdW functional was em¬ 
ployed Hanoi- This vdW density functional has been 
applied to a wide range of materials and proven to be of 
good accuracy [27H33] . The Brillouin zone sampling was 
performed on a T-centered 3x3x1 k-grid, with a smear¬ 
ing of 0.1 eV using the method of Methfessel and Pax¬ 
ton to first order [34]. To ensure good accuracy for the 
Bader-partitioning scheme, the electronic charge density 
was calculated on a dense mesh of 432 x 432 x 448 points 
in the simulation cell. Electronic energies were converged 
to 10“^ eV and the ionic relaxations were stopped after 
converging forces between the interacting atoms to better 
than 0.01 eV/A. Eollowing the investigation by Garhofer, 
who also provided us with initial structural data m , we 
choose a plane wave cutoff of 300 eV, which is the min¬ 
imum recommended value for carbon, and at the same 
time larger than the suggested value for iridium and tung¬ 
sten. 

To partition the electronic charge density p in our sim- 
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(a) Side view 



(b) Top view 

FIG. 1: Side (a) and top (b) view of a tungsten tip on 
graphene/Ir(lll). Iridium atoms are shown in yellow, 
carbon in brown and tungsten in grey. 


ulation cell into single atoms we use Bader’s quantum 
theory of atoms in molecules (QTAIM) [35]. In contrast 
to other similar approaches [36H4T] . Bader’s method pro¬ 
duces non-overlapping atomic domains with well-defined 
boundaries, which are perfectly suited to analyze the con¬ 
tact between two adjacent bodies. The necessary and suf¬ 
ficient condition that needs to be fulfilled to define the 
boundaries of a selected atom according to the QTAIM 
is formulated using the basic quantity in DFT, namely 
the electronic charge density p and is given as 02], 

Vp(r) • n(r) = 0 Vr G S{y) . (1) 

Here S is the boundary surface of the atom and n(r) is 
the unit vector normal to this surface. The condition 
states that the flux of the gradient field of the charge 
density, Vp, through the boundary surface S must van¬ 
ish, which is thus called a zero flux surface. The Bader 
analysis in this work was performed with the code devel¬ 
oped by Henkelman, Sanville, and Tang which is directly 
compatible with the format of VASP-output files [43ll45] . 


RESULTS 

If one seeks to define the real contact area on an atomic 
scale, it is quite natural to think about the size, shape, 
and deformations of the involved atoms. Once the size 
and shape of all atoms in the contact region are deter¬ 
mined, the calculation of the real contact area is reduced 


to a simple summation of the regions that are in contact, 
provided one can distinguish unambiguously between the 
two contacting bodies. Since p formally is non-zero every¬ 
where, the Bader-atoms at the surfaces of the contacting 
bodies extend into the vacuum region to infinity or un¬ 
til they encounter another atom. This would mean that 
contact between two bodies is established at all distances, 
which is a clearly unphysical result, unless one defines a 
density cutoff. This density cutoff pcut cannot be chosen 
arbitrarily, since it directly influences the contact area. 

A possibility to extract a value for pcut is to analyze 
the interaction potential between the tip and the sur¬ 
face, divide it into a long- and a short-range part and 
define the contact at the onset of the short-range in¬ 
teraction, analogous to Mo, Turner, and Szlufarska [10]. 
This procedure defines a density cutoff pcut so that the 
Bader-partitioning yields contact only after the onset of 
short-range interactions. However, as long-range interac¬ 
tions are included implicitly in the exchange-correlation 
potential that we use (see section Computational Meth¬ 
ods), the separation into a long- and a short-range part is 
not straightforward. A more unambiguous way to define 
the onset of contact is to analyze the atomic relaxations 
that happen if the tip is lowered towards the surface. We 
distinguish the distance for the static, unrelaxed system 
{ds) and the relaxed distance (d^), which are both mea¬ 
sured between the tip apex atom and the carbon atom 
directly beneath it. For large distances no relaxations 
will happen, although there might be attraction due to 
vdW forces, and the distance dr in the relaxed system 
will be equal to the (static) distance ds measured before 
relaxing the system. At some point during the approach 
of the tip stronger forces will cause relaxations, which 
will result either in a “snap” or “jump” to contact (a 
phenomenon often observed in AFM experiments) if the 
interaction is attractive, or in a depression of the surface 
layer if the interaction is purely repulsive. A sketch of 
this process is given in figure [^ In any case, the dr ver¬ 
sus ds curve will have a discontinuity at some distinct 
distance where the system begins to strongly interact. 
Below this distance the system will try to hold the ideal 
distance between tip and surface. It is straightforward 
to identify this discontinuity as the onset of contact. 

In the examined system the interaction between the 
tip and the surface is attractive at the onset of contact 
and a jump to contact occurs between ds = 3.65 A and 
ds = 3.53 A (see figure 3a), such that dr is changing 
from 3.5 A to 2.7 A, accordingly. Most of the movement 
is done by the surface, which jumps upwards to meet 
the tip (see figure 3b). This means that when the tip 
support is lowered by only 0.12 A, the distance between 
the tip apex and the surface is reduced by 0.8 A. This 
allows us to define the onset of contact at 3.6 A and 
to tune the value of the density cutoff p^t accordingly. 
Note that an infinitely stiff cantilever is assumed in our 
calculations, as the uppermost atoms of our tip are kept 
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FIG. 2: Sketch of possible tip-surface interactions. If 
the tip is far away from the surface (middle panel), 
atomic relaxations will not have an effect on the 
distance between the tip and the surface. If the tip gets 
closer, the surface will interact with the tip and will 
either jump towards the tip, if the interaction is 
attractive (left panel), or will get depressed, if the 
interaction is repulsive (right panel). This effect can be 
used to define the onset of contact. 


rigid. In principle, the influence of a more compliant 
AFM apparatus on the initial jump to contact, however, 
could also be modeled by using a multiscale approach. 

Once pcut is selected the determination of the real con¬ 
tact area for each distance is straightforward. Since the 
partitioning code produces only Bader volumes rather 
than zero flux surfaces, we have to construct the contact 
area from these data. To this end the Bader-volumes 
of both contacting bodies are added up and by pairwise 
comparison of the respective values for neighboring grid 
points a point cloud forming the contact area is gener¬ 
ated. The contact area can now be obtained by triangu¬ 
lation. 

We calculated contact areas for cutoff densities from 
10“^e/A^, which is the default cutoff in the partition¬ 
ing code by Henkelman, Sanville and Tang [43II45] . up 
to a cutoff of 10“^e/A^. While the default value of 
pent = 10“^ e/A? and other low cutoff densities are giving 
sizable contact areas for all distances, we approach the 
desired effect of establishing contact only for dg < 3.6 A, 
for a value of peut r\j 10 ^e/A^. Obviously, a very high 
Peut is unphysical, as the number of electrons that are 
“lost” into the vacuum region increases with rising cut¬ 
off. This means that we want to select a value that is 
high enough to guarantee that the contact area A is only 
non-zero after the snap to contact has occurred, but is 
otherwise as low as possible. We analyzed several values 
of peut ranging from 7.5 x 10“^ e/A^ to 1.0 x 10“^ e/A^ in 
order to find the lowest value that still satisfies these con¬ 
ditions, resulting in an optimal value of peut = 3 x 10“^ 
electrons per A^. 

Of course, changing the contact site of the tip away 
from a position directly above a carbon atom or to an¬ 
other section of the moire structure could conceivably 
change the exact point of the jump to contact and thus 
modify the value of Peut- However, the obtained cutoff 
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FIG. 3: (a) Distance dr in the relaxed system (red 
crosses) versus the (static) distance dg between a rigid 
tungsten tip approaching graphene/Ir(lll). The dashed 
line gives dr = dg. Attractive interactions cause a jump 
to contact between dg = 3.65 A and dg = 3.53 A, which 
is marked with a vertical dotted line. After the jump to 
contact the relation between dr and dg is also linear 
(solid black line). This line crosses the dr = dg line at 
the equilibrium point, where the graphene layer has 
relaxed back into its original shape. A hypothetical 
curve for a purely repulsive interaction is sketched by 
the dashed-dotted line to illustrate that the method is 
also viable if no jump to contact is occurring in the 
system, (b) Displacement Dc of the carbon atom 
situated directly below the tip apex with respect to its 
initial position, versus the distance dg. The jump to 
contact is again marked by a vertical dotted line. 


density of peut = 5 x 10“^ is low enough to result in 
an essentially flat graphene surface and hence changing 
the tip position should only marginally alter the com¬ 
puted contact area. As this paper is mainly concerned 
with the presentation of a new approach in defining and 
calculating the real contact area ah-initio^ and given the 
rather large computational effort [46] , we decided against 
repeating our calculations on different contact sites and 
calculating an average. 

In a preliminary calculation of the same tip on an fee 
copper (111) surface we found an optimal charge den¬ 
sity cutoff value of 5.3 x 10“^ electrons per A^, which 
is approximately the same as for the graphene/Ir(lll) 
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system. 

As the optimized cutoff of 5 x 10“^e/A^ is 50 times 
larger than the default value it is important to check if it 
is still a reasonable number and does not give any unphys¬ 
ical results. We therefore calculated the nominal number 
of electrons that are assigned to the vacuum region and 
thus are not part of any Bader-atom. For the default pcut 
of 1 X 10 ^ e/A^ only about half of an electron is not rep¬ 
resented by a Bader-atom. For the cutoff value needed 
for the calculation of the contact area, 5 x 10 ^e/A^, 
this number is increased to nearly 30 electrons. Although 
this value seems to be very large, one has to consider the 
total system size, which includes 3776 electrons. Thus, 
the relative number of “missing” electrons is below 0.8%. 
We also evaluated the Bader-radii Rb of a single tung¬ 
sten atom and a single carbon atom in a box. The value 
for tungsten of 2.9 A obtained for pcut = 5 x 10“^ elec¬ 
trons per A^ is more than twice as large as the empirical 
atomic radius, 1.35 A m, and about 1 A larger than the 
calculated atomic radius of 1.93 A [48]. Reported values 
for the atomic radius of carbon reach from the calculated 
value of 0.67 A [48], over the empirical value of 0.70 A [47] 
to a van der Waals radius of 1.70 A [49]. As for tungsten, 
also the carbon Bader-atom radius for a cutoff density of 
5 X 10“^ electrons per A^ is significantly larger than all 
these reported values at 2.30 A. This ensures that neither 
the tip atoms nor the surface Bader-atoms are artificially 
small. It also shows that it is questionable to assume that 
contact between two bodies is established only after the 
surface atoms overlap if one uses spherical atoms and 
traditional radii. 

It is worthwhile to compare our approach to the onset 
of contact with the method by Mo, Turner, and Szlu- 
farska m, which uses the beginning of short-range inter¬ 
action as a criterion for contact. As already mentioned, 
the distinction between long-range and short-range forces 
is not trivial, but it might be approximated by disabling 
the long-range contributions in the correlation potential 
of the optB86b-vdW functional. We calculated the cor¬ 
responding energies at the vdW relaxed positions and 
fitted the data with a Morse function m- The interac¬ 
tion strength of this short-range potential at the jump to 
contact is 2.0% of the total potential depth which could 
be classified as the “beginning of the interaction”. This 
means that our approach, at least for the system investi¬ 
gated here, is in accordance with the approach of refer¬ 
ence uni However, an interaction strength of 1%, 5% or 
even 10% of the short-range binding energy could also be 
reasonably selected as the “beginning of the interaction”, 
each leading to different results. This highlights the ad¬ 
vantages of using the jump to contact as the criterion for 
the initial point of contact, as no further assumptions are 
needed proceeding in this manner. 

Figure shows a decomposition of the contact area 
into contributions of the tip apex (one atom) and con¬ 
tributions from the second tip layer (four atoms). We 



relaxed distance d^. [A] 


FIG. 4: Ah-initio real contact area A obtained for 
Pcut = 5 X 10“^ e/A^ versus the distance dr in the 
relaxed system between the tungsten tip and 
graphene/Ir(lll). Blue crosses give the contribution of 
the tip from the apex atom and green plus signs show 
the contribution from the second layer (four atoms). 
The total contact area (red circles) is given by the sum 
of these two contributions and is fitted by an 
exponential (dashed line). The inset shows the total 
contact area versus the static distance dg. The solid line 
is an exponential function resulting from equation ([^ 
and the linear relation between dr and dg (see figure [3^ . 
The dotted vertical line marks the jump to contact. 


find that the second layer only contributes to the total 
contact area for the three closest distances but is then re¬ 
sponsible for nearly all of the increase. The dashed line 
in figure is an exponential fit of the form 

A{dr) = , ( 2 ) 


to the 7 non-zero data points (red circles) with the co¬ 
efficients 4.2 A“^ and 3.0 A. The factor 

Aa = lA^ is included for dimensional reasons and has 
not been used as a fitting parameter. Although the ex¬ 
ponential fit is not perfect, the agreement with the data 
is certainly reasonable, especially considering that only 
two fitting parameters were used. Also the point of van¬ 
ishing contact is predicted well, although only points of 
positive contact area were considered for the fit. As there 
is a linear relation between dr and dg in the region where 
contact is established {dr = i^dg 6 = 0.23dg -h 1.84; 
see solid black line in figure 3a), it is also possible to 
express the contact area A through the static distance 
dg , which is easier accessible in experiments through the 
vertical displacement of the tip support. In the relation 
A{dg) = Aa exp [—As (dg — Ag)], the decay constant is 
smaller than in A{dr) with A^ = XrK. 1.0 A“^, while 
As = (A^ — 5)/n 5.0 A is increased compared to A^, 
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and Aa = 1 is the same dimensionality factor as be¬ 
fore. This relation is plotted in the inset of figure 
Quite naturally only the region after the jump to contact 
is represented well. 

Figure shows the geometrical shape of the non¬ 
vanishing real contact area for four distances. The differ- 



the iridium substrate, the contact area is beginning to 
show a pronounced bowl shape (figure [^, which in¬ 
creases in depth for decreased distance (figure [5d|). Note 
that for the closest distance (figure [5d| ds = 1.30 A) not 
only the threefold symmetry of the graphene layer is vis¬ 
ible in the center, but the edges of the bowl have the 
fourfold symmetry of the second tip layer. 

We can also analyze how many carbon atoms are in 
contact with the tip for each distance and compare the 
contact area predicted by our approach with the results 
by Mo, Turner, and Szlufarska m- To this end we count 
every surface Bader atom that touches our tip as contact¬ 
ing, a different approach than described in the introduc¬ 
tion, since we have no pairwise forces at our disposal. 
In our case, with a cross section of the simulation cell 
Ac r\j 262 and 200 carbon atoms in the graphene 
layer the contribution per atom to the real contact area 
is Aat = Ac/200 = 1.31 A^. Each carbon atom has 3 
nearest neighbors in dnn = 1.42 A, 6 next nearest neigh¬ 
bors at 2.46 A, 3 third nearest neighbors at 2.84 A, and 
6 fourth nearest neighbors at 3.76 A distance. Already 
directly after the jump to contact at dg = 3.53 A, more 
than one carbon atom is in contact with the tip, although 
the majority of the contact is formed by the central car¬ 
bon atom which is responsible for 5.90 A^ of the total 
6.82 A^. This is about 30% more than the contact area 
predicted in reference [TOl with 4Aat = 5.25 A^. The area 
contributed by the central carbon atom alone exceeds 
the value of 4Aat by ^ 12%. The next nearest and third 
nearest neighbors begin to play a role at dg = 2.18 A, 
contributing to about 7% of the total area of 16.89 A^. 
Here the method by Mo, Turner, and Szlufarska gives 
a very comparable area of 13Aat = 17.05 A^. However, 
the 9 outermost atoms that contribute 7% to the con¬ 
tact area in our approach are responsible for nearly 70% 
of the contact area in reference [10] considering all con¬ 
tacting atoms equally. The situation at dg = 1.85 A is 
visualized in figure with the method of reference m 
still giving an area of 13Aat = 17.05 A^, while our ap¬ 
proach yields 21.13 A^. Only for the two closest positions 
at dg = 1.51 A and dg = 1.30 A more than 13 atoms are 
in contact, according to the Bader partitioning, and the 
central 13 are still responsible for 98% and 87% of the 
contact area, respectively. Including the 6 fourth near¬ 
est neighbors into the model by Mo, Turner, and Szlu¬ 


farska [To], leads to 19Aat = 24.92 A^ for both of this 
distances, while our approach gives A(1.51) = 28.89 A^ 
and A(1.30) = 35.03 A^. Thus, the results are compara¬ 
ble, but the outermost atoms are again over represented 
compared to our approach. Our model offers higher res¬ 
olution of the real contact area and allows for a distance 
dependent contribution of each atom. It is important to 
note that our contact areas are curved and have a more or 
less pronounced bowl shape while Mo, Turner, and Szlu¬ 
farska consider flat contact areas (see figure]^. Overall 
both methods show fair agreement. 

Our chosen system, which has been proven to accu¬ 
rately model the interaction between a tungsten tip and 
moire graphene on Ir(lll) [T7|, limits our investigation 
to the attractive region {dr > 2.24 A) and small positive 
loads (2.18 A < dr < 2.24 A). For dr < 2.18 A, the tip 
forms bonds with the iridium substrate leading again to 
negative values of the load. Thus, it is difficult to predict 
the behavior of the real contact area A dependent on the 
load L. However, we can assume that the interaction po¬ 
tential E{dr) can be approximated by a Morse potential 
in the vicinity of the minimum m, 


Em {dr) = Eq 


{ 


I _ ^-lidr-do) 



(3) 


where Eq = 2.33 eV is the depth of the potential at the 
equilibrium position do = 2.24 A, which we can get di¬ 
rectly from our data. Thus only 7 has to be fitted, re¬ 
sulting in 7 = 4.11 A“^. We can now derive the load 
L = —dEM/ddr yielding 


L {dr) — —2jEo 


I _ ^-jidr-do) 


^-^(dr-do) 


(4) 


Solving this for dr produces 

dr = do - , ( 5 ) 
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where the dimensionless function ^{L) is 

/-/T\ 1 E y/l — 4:U . ^ L 

^{L) = -^- , with u=— . (6) 

Now it is possible to express the real contact area depen¬ 
dent on load using equations (© and (§, which provides 
a power law, 

. ( 7 ) 

As 7 = 4.11 A“^ and = 4.19 A“^ the exponent is very 
close to 1, thus we arrive at a linear dependence of A on 
^(L), namely A{L) = C^{L) with C = 24.15 A^, which 

corresponds to an increase of A with L to the power of 

1 

2 • 

While we believe that our ab-initio approach using the 
QTAIM for calculating the real contact area is intuitive 
and accurate, its limitations have to be discussed as well. 
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(c) ds = 1.85 A 



x[A] 

(d) ds = 1.30 A 


FIG. 5; Contact formed by lowering a tungsten tip onto a graphene/Ir(lll) surface. The contact area increases from 
(a) 6.8 A2 , over (b) 9.1 A^, and (c) 21.13 A^, to (d) 35.0 A^. Different depths of the curved contact areas are coded by 
color contours with the lowest value set to zero. Please note the different color bars and axes scaling in each panel. 


As the charge density is required to perform the Bader 
partitioning our approach is limited to system sizes where 
ah-initio calculations are still feasible. This limits us to a 
single asperity case at the moment, where only a handful 
of atoms interact with the surface. However, since our 
system was used to explain some experimental results m 
[m, we are confident that our scheme is applicable to 
real systems, albeit only for sharp tips and low loads. 
With the ever increasing power of modern computers and 
better scaling codes, on the other hand, it might soon be 
feasible to calculate systems with thousands or millions 
of atoms and hence to study interesting phenomena such 
as multi-asperity contacts. Thus, our approach might be 
used in MD calculations as well to directly investigate the 
influence of the real contact area on frictional forces m- 


CONCLUSION 

We propose a new ab-initio approach for calculating 
the real contact area between a tip and a surface. We 


apply Bader’s quantum theory of atoms in molecules 
(QTAIM) to determine the volumes and shapes of the 
atoms in contact together with their contact areas at 
each given distance. We define a specific density cutoff 
Pcut foi* this partitioning to confine the Bader volumes to 
realistic values. This cutoff density is obtained by using 
the discontinuity in the dr versus dg curve to define the 
initial point of contact in a perspicuous way, which in 
the examined system occurs due to a jump to contact, 
commonly observed in AFM experiments. This defines a 
lower bound for the cutoff density which is then the opti¬ 
mal value, since pcut needs to be minimized to include the 
maximum number of electrons. Thus, our approach re¬ 
mains essentially ab-initio^ as the only parameter needed 
can be determined from properties of the system. We 
believe that the jump to contact is a less ambiguous way 
to define the onset of contact than using a partitioning 
of the interaction in long- and short-range regions HO], 
or equating contact with repulsive interactions min]. 

For decreasing the real tip-sample distance dr an ex¬ 
ponential increase of the real contact area A is found. 































FIG. 6: Top view (left) and side view (right) of the real 
contact area resulting from our ab-initio approach using 
Bader atoms (color code depend ing on height) for a 
distance of dg = 1.85 A (see figure 5c), compared to the 
flat contact area from the model by Mo, Turner, and 
Szlufarska (green) [TO], Carbon and tungsten atoms are 
sketched as red and black dots, respectively. 


This is a combined effect of the jump to contact and the 
preferred distance of the tip apex relative to the surface 
atom below it, which first jumps up to meet the tip and 
then gets pressed below its equilibrium position for closer 
separations. As dg is linear dependent on dr^ we can also 
express the exponential relation A{dr) through which 
can be better controlled in experiments than d^. 
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